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Abstract 

In this short note, by analyzing the dispersion and dissipation of explicit finite 
difference scheme, a dispersion-dissipation condition is derived to determine the 
minimum dissipation required to damp the artificial high-wavenumber waves in 
the solution. The example application to our previous developed WENO-CU6-M2 
scheme suggests that this condition can be used as an general guidance on optimizing 
the dissipation of a numerical scheme. 
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1 Introduction 

The spectral analysis of finite difference scheme shows that the dispersion and 
dissipation introduced by the scheme determine the errors related to the prop- 
agation of a simple wave with given wavenumber. While the original spectral 
analysis is only for linear schemes, Pirozzoli 9| developed a numerical spec- 
tral analysis for general non-linear finite difference schemes, such as shock- 
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ll| found that the dispersion and 



capturing scheme. Recently, Sun et al. 
dissipation of a class of the explicit linear schemes can be controlled inde- 
pendently. Since the pioneer works of Tarn and Webb Jjj] on computational 
aero-acoustics the spectral analysis has been applied widely to the design of 
high-order numerical schemes, especially those for the direct numerical sim- 
ulation (DNS) of turbulent flow, for accurately resolving small scale in both 



amplitude and phase 



15 



10 



labj. 



While it is generally accepted that the dispersion of a numerical scheme should 
be minimized according to some chosen criteria, there are no general guidelines 
on how the dissipation should be optimized [11 1 . For example, it is known 
that the weighted essentially non-oscillatory (WENO) schemes using upwind 
schemes as the underlying optimal linear scheme are too dissipative for DNS. 
However, the minimum dissipation produced by the central scheme may not 
be sufficient to suppress numerical oscillations 6|. Pirozzoli |8| pointed out 
that a certain amount of dissipation is desirable to damp as much as possible 
the high-wavenumber waves with incorrect propagation speed. The difficulty 
is to determine how much dissipation a numerical scheme should provided to 
damp the artificial high-wavenumber waves, but not the physical waves to be 
resolved. 



The objective of this short note is to show that the minimum sufficient dis- 
sipation of a numerical scheme is actually determined by its dispersion with 
a dispersion-dissipation condition. This condition is then applied to alleviate 
the difficult of artificial waves produced by our previous developed WENO- 
CU6-M2 scheme [2j. Two straightforward modifications are proposed together 
with the numerical example on the two-blast wave interaction problem. 
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2 Dispersion-dissipation condition 



For completeness, we begin from the spectral analysis of the one- dimensional 



linear advection equation 



13| 



du du q 
dt dx 



-oo < x < +00, u(x, 0) = e v 



with the initial condition of a monochromatic sinusoidal with wavenumber w 
and unit amplitude where c is the advection speed. Eq. ([!]) is discretized in 
the spatial domain such that Xj = jAx, j = 0, N, where Ax is spacing and 
Uj = u(xj) is the node value, and the semi-discretized form by the method of 
lines yields a system of ordinary differential equations 



0. 



N, Uj{0) = e^*, 



(2) 



where 6 = wAx is the reduced wavenumber, and w' ~ 



du 

3 ~ dx 



is an finite 



difference approximation. If an explicit linear scheme is applied 



l=-q 



(3) 



the analytical solution of Eq. (j2J) is 



(4) 



where 



r r 

H<f>) =-»E a ' e ^ =E«i t sin (^) - icos^)] , 



(5) 



is the modified wavenumber. If Eq. (j3J) is compared with the solution of the 
exact Eq. (CQ) at the Xj 

Uj (t) = e tw ^ Ax ~ ct \ (6) 

it is easy to find that the introduction of the finite difference approximation 
Eq. 02]) leads two type of numerical errors. The first type of errors, called 
dissipation, modifies the amplitude of the wave with the factor e~ xt , where 
A = iIm($(0))c/Ax is the decay rate. In the classical theory of linear stability, 
if a > 1 or A < 0, the dissipative error leads to the exponentially growth of 
the amplitude and the semi-discretization Eq. ([2]) is unstable. The second 
type errors, called dispersion, modifies the propagation speed of the wave to 
cc*((j)), where c*(<f>) = Re($(0))/0 is the normalized phase velocity. Usually, 
for a complex wave packet with a wavenumber spectrum, the apparent velocity 
of a wave component with phase velocity c*(4>) is given by the group velocity 
V = cd(4>c* (0)) / d(p = cRe($(0))', Note that, while the wave is spreading 
away from the main wavefront with the speed |Ac| = \V — c\, its magnitude 
is decreasing with the decay rate A. Therefore, the distance 5 from the main 
wavefront in which the wave can survive is 

6 ~ A ^Im($(l))c- (7) 

If the solution is free of artificial waves, i.e., S ~ Ax, the condition must be 
satisfied is the dispersion-dissipation ratio 

|Re W) y-l| 

ilm($(0)) ' 1 ; 

if the dispersion error is not negligible. Since 5 decreases with increase of A, 
Eq. OH]) actually determines the minimum sufficient dissipation required to 
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damp the artificial waves. Specifically, for the explicit linear scheme Eg. (13"]) 
the above condition is 



EL-o a ll COS(Z0) - 1 



E/=-fl a l COS(Z0) 



1. 



(9) 



Note that, for non-linear schemes, since the modified wavenumber Eq. §5§ 



can be obtainec 
relation (ADR) 
Eq. (JED- 



approximated numerically by the approximated dispersion 
9j, one can still define the dispersion dissipation condition by 



3 Application to WENO-CU6-M2 scheme 

The WENO-CU6-M2 scheme 2J is an adaptive central-upwind weighted es- 
sentially non-oscillatory (WENO) scheme with 6th-order accuracy for smooth 
solutions. The scheme uses a 6-point stencil to achieve adaptively, by a non- 
linear weighting procedure, either the approximation of the optimal 6th-order 
central scheme in smooth regions of the solution, or a numerical stable 3rd- 
order approximation when there are discontinuities or shocks in the solution. 
The nonlinear weights are given by 

" r = T£j? °' = * P + + ^ ft, r + xA^ j ' (10) 

where d r = {0.05, 0.45, 0.45, 0.05} are the optimal weights leading to the ap- 
proximation of the 6th-order central scheme, q = 4 is an integer parameter, 
C q = 10 3 is a positive constant parameter and e = 1/x = 10~ 8 is a small 
positive number. For the WENO methodology and the details of computing 
p3,n P3,ave and tq, please refer to Refs. 4j, y]. As shown in Fig. [T^,, the mod- 
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ified wavenumber of WENO-CU6-M2 is very close to that of the 6th-order 
central scheme until about (f) = tt/2. Since the latter produces considerable 
dispersion beyond 0=1 but without dissipation, as shown in Fig. [Tb, the 
dispersion-dissipation ratio of WENO-CU6-M2 does not satisfy the condition 
Eq. (jSJ) in the wavenumber range aroung <ft = 1.1 to n/2, and has the max- 
imum value about r = 100 at <p = 1.302. Therefore, WENO-CU6-M2 may 
produces artificial waves in such wavenumber range. 

To illustrate the artificial waves in the solution produced by WENO-CU6-M2, 



the two-b 
Colella |l| 



ast-wave interaction problem, which is taken from Woodward and 
, is simulated on 400 grid points. The reference "exact" solution 



is a high-resolution numerical solution on 3200 grid points calculated by the 
WENO-CU6 scheme 3j. The numerical solution in Fig. [2] shows that, while 
WENO-CU6-M2 achieves very high resolvability for the density valley at about 
x = 0.75 and peak at about x = 0.78 (see Fig. [2^ and b), small artificial waves 
can be observed in the close-up views of the velocity profile (see Fig. [2b and 
d). These waves have wavelengths about 4 to 6 grid points, which corresponds 
to the wavenumber range with large r as shown in Fig. [Tb. 

The difficulty of artificial waves can be alleviated either by increasing the dis- 
sipation (denoted as modification A), or by modifying both dispersion and dis- 
sipation (denoted as modification B), to satisfy the dispersion-dissipation re- 
lation of Eq. (jSJ). For modification A, other than the 6th-order central scheme, 
the 5th-order linear scheme with the constraint r = 10 at <fi = 1.302 according 
to Eq. is chosen as the optimal scheme. The resulting the optimal weights 
are d r = {0.065, 0.495, 0.405, 0.035}. As shown in Fig.UJi and b, such modifica- 
tion dos not change the dispersion, but achieves r < 10 in the entire wavenum- 
ber range. The numerical solution of the two-blast-wave interaction problem, 
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as shown in Fig. [2^ and b, suggests no notable differences from those obtained 
by the original scheme. However, the artificial waves in the velocity profile are 
eliminated (see Fig. and d). For modification B, a 4th-order scheme, whose 
dispersion is optimized by minimizing the integral error E = Jq e 6 ^ n ~^\c*\ 2 d(f) 
with the constraint c* < 0.015 for all <fi < <f) c , where C = 1.5 is a critical 
wavenumber, and whose dissipation is separately constrained by r = 10 at 
= 7r/2 111 ], is chosen as the underlying optimal scheme. The resulting opti- 
mal weights are d r = {0.0904546,0.4440908,0.3922727,0.0731819}. As shown 
in Fig. [T^,, the modified scheme achieves both considerable improvement on 
dispersion, and r < 10 in the entire wavenumber range. Again, as shown in 
Fig. [2] the numerical solution obtained by the modified scheme suggests no 
notable differences from those by the original scheme, and the artificial waves 
in the velocity profile are eliminated. 



4 Concluding remarks 



In this short note, we have derived the dispersion-dissipation condition which 
determines the minimum dissipation that a numerical sufficient finite differ- 
ence scheme required to damp the artificial high-wavenumber waves in the 
solution. With the help of approximated dispersion relation, the condition can 
also be defined for non-linear schemes. The example application to WENO- 
CU6-M2 scheme suggests that this condition can be used as an general guid- 
ance on optimizing the dissipation of a numerical scheme. 
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Figure 1. Approximated dispersion relation of various numerical schemes: (a) modi- 
fied wavenumber and (b) dispersion-dissipation ratio computed with r = ^m^fi))+e , 
where e = 10~ 3 . 
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Figure 2. Numerical solution of the two-blast- wave interaction problem at t = 0.038 
on a 400 points grid: (a) density and (b) velocity profiles, (c) and (d) close-up views 
of the velocity profile. 
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